#define _SIZE_T_DEFINED
#ifndef __CUDACC__
#define __CUDACC__
#endif
#ifndef __cplusplus
#define __cplusplus
#endif

#include <cuda.h>
#include <device_launch_parameters.h>
#include <texture_fetch_functions.h>
#include "float.h"
#include <builtin_types.h>
#include <vector_functions.h>
static __inline__ __host__ __device__ float4 fminf(float4 a, float4 b)
{
	return make_float4(fminf(a.x,b.x), fminf(a.y,b.y), fminf(a.z,b.z), fminf(a.w,b.w));
}
static __inline__ __host__ __device__ float4 fmaxf(float4 a, float4 b)
{
	return make_float4(fmaxf(a.x,b.x), fmaxf(a.y,b.y), fmaxf(a.z,b.z), fmaxf(a.w,b.w));
}

extern "C"  
{
	typedef struct {
		float x, y, z;
	} Point;

	typedef struct {
		float x, y, z;
	} Vector;

	typedef struct {
		Point o;
		Vector d;
		float mint, maxt;
	} Ray;

	typedef struct {
		float Distance;
		float U, V; // Barycentric coordinates of the hit point
		unsigned int Index;
	} RayHit;

	typedef struct {
		unsigned int v[3];
	} Triangle;

	typedef struct {
		Point pMin, pMax;
	} BBox;

	typedef struct {
		BBox bbox;
		unsigned int primitive;
		unsigned int skipIndex;
	} GpuBVHNode;

	typedef struct {
		BBox bbox;
		int primOffset; //+4 
		int secondChildOffset;
		unsigned char nPrimitives; //+1
		unsigned char axis; //+1
	} GpuLBVHNode;

	__device__ float4 cross(const float4 &a, const float4 &b) {
		return make_float4
			(a.y * b.z - a.z * b.y,
			a.z * b.x - a.x * b.z,
			a.x * b.y - a.y * b.x,
			0.f);
	}

	__device__ float dot(const float4 &a, const float4 &b)
	{
		return a.x*b.x+ a.y*b.y+a.z*b.z +a.w*b.w;
	}

	__device__ float4 operator+(const float4 &a, const float4 &b) {

		return make_float4(a.x+b.x, a.y+b.y, a.z+b.z, a.w+b.w);
	}
	__device__ float4 operator*(const float4 &a, const float4 &b) {

		return make_float4(a.x*b.x, a.y*b.y, a.z*b.z, a.w*b.w);
	}
	__device__ float4 operator/(const float &a, const float4 &b) {

		return make_float4(a / b.x, a / b.y, a / b.z, 0.f);
	}
	__device__ float4 operator-(const float4 &a, const float4 &b) {

		return make_float4(a.x-b.x,a.y-b.y,a.z-b.z,a.w-b.w);
	}

	__device__ void TriangleIntersect
		(
		const float4 rayOrig,
		const float4 rayDir,
		const float minT,
		float *maxT,
		unsigned int *hitIndex,
		float *hitB1,
		float *hitB2,
		const unsigned int currentIndex,
		const Point *verts,
		const Triangle *tris
		)
	{
		const Point *p0 = &verts[tris[currentIndex].v[0]];
		const Point *p1 = &verts[tris[currentIndex].v[1]];
		const Point *p2 = &verts[tris[currentIndex].v[2]];


		float4 v0 = make_float4(p0->x, p0->y, p0->z,0.f);
		float4 v1 = make_float4(p1->x, p1->y, p1->z,0.f);
		float4 v2 = make_float4(p2->x, p2->y, p2->z,0.f);

		float4 e1 = v1 - v0;
		float4 e2 = v2 - v0;
		float4 s1 = cross(rayDir, e2);

		const float divisor = dot(s1, e1);
		if (divisor <= 0.f)
			return;
		const float invDivisor = 1.f / divisor;

		// Compute first barycentric coordinate
		float4 d = rayOrig - v0;
		float b1 = dot(d, s1) * invDivisor;
		if (b1 < 0.f) return;

		// Compute second barycentric coordinate
		const float4 s2 = cross(d, e1);
		const float b2 = dot(rayDir, s2) * invDivisor;
		if (b2 < 0.f) return;

		const float b0 = 1.f - b1 - b2;
		if (b0 < 0.f)  return;

		// Compute _t_ to intersection point
		const float t = dot(e2, s2) * invDivisor;
		if (t < minT || t > *maxT)  return;
		*maxT = t;
		*hitB1 = b1;
		*hitB2 = b2;
		*hitIndex = currentIndex;
	};
	__device__ int BBoxIntersectP(
		const float4 rayOrig, 
		const float4 invRayDir,
		const float mint, const float maxt,
		const float4 pMin, const float4 pMax) 
	{
		const float4 l1 = (pMin - rayOrig) * invRayDir;
		const float4 l2 = (pMax - rayOrig) * invRayDir;
		const float4 tNear = fminf(l1, l2);
		const float4 tFar =  fmaxf(l1, l2);

		float t0 = max(max(max(tNear.x, tNear.y), max(tNear.x, tNear.z)), mint);
		float t1 = min(min(min(tFar.x, tFar.y), min(tFar.x, tFar.z)), maxt);

		return (t1 > t0);
	}
	__device__ static inline int IsDirNegative(const float4 dir, const unsigned char axis)
	{
		if (axis == 0)
			return dir.x < 0;
		if (axis == 1)
			return dir.y < 0;
		if (axis == 2)
			return dir.z < 0;
		return 0;
	}

	__global__ void IntersectLBvh(
		Ray* rays,
		RayHit *rayHits,
		Point *verts,
		Triangle *tris,
		GpuLBVHNode *bvhTree,
		const int triangleCount,
		const int nodeCount,
		const int rayCount
		)
	{
		const int gid = blockIdx.x*blockDim.x+threadIdx.x;//blockDim.x * blockIdx.x + threadIdx.x
		if (gid >= rayCount)
			return;


		float4 rayOrig,rayDir;

		float minT, maxT;
		{
			float4 *basePtr =(float4 *)&rays[gid];
			float4 data0 = (*basePtr++);
			float4 data1 = (*basePtr);

			rayOrig = make_float4(data0.x, data0.y, data0.z, 0.f);
			rayDir =  make_float4(data0.w, data1.x, data1.y, 0.f);

			minT = data1.z;
			maxT = data1.w;
		}
		
		float4 invRayDir =  make_float4(1.f / rayDir.x , 1.f / rayDir.y ,1.f  / rayDir.z ,1.f);

		unsigned int todoOffset  = 0;
		unsigned int hitIndex = 0xffffffffu;
		unsigned int currentNode = 0; // Root Node
		float b1, b2;
		int todo[32];
		//const unsigned int lid = get_local_id(0);
		todo[0] = 0;
		float4 pMin, pMax, data0, data1;
		float4 *basePtr;
		unsigned int indexOfHit = 0xffffffffu; // Non-existent
		float dist = 1e17f, u, v;
		while (true) {
			basePtr =(float4 *)&bvhTree[currentNode];
			data0 = (*basePtr++);
			data1 = (*basePtr);

			pMin = make_float4(data0.x, data0.y, data0.z, 0.f);
			pMax = make_float4(data0.w, data1.x, data1.y, 0.f);

			if (BBoxIntersectP(rayOrig, invRayDir, minT, maxT, pMin, pMax)) 
			{
				
				const int triangleIndex = bvhTree[currentNode].primOffset;
				const char count = bvhTree[currentNode].nPrimitives;

				//indexOfHit = triangleIndex;
				//dist = maxT;
				//break;

				if (count > 0)
				{
					for(int ti=0;ti<count;ti++)
					{
							TriangleIntersect(rayOrig, rayDir, minT, &maxT, &hitIndex, &b1, &b2, triangleIndex, verts, tris);
							if (maxT < dist) 
							{
								rayHits[gid].Distance   = maxT;
								rayHits[gid].U = b1;
								rayHits[gid].V = b2;
								rayHits[gid].Index = hitIndex;
							}
					} 
					if( todoOffset == 0 ) {
						break;
					}
					currentNode = todo[--todoOffset];
				} 
				else 
				{
					const unsigned char axis = bvhTree[currentNode].axis;
					if( IsDirNegative(invRayDir, axis) ) 
					{
						todo[todoOffset++] = currentNode + 1;
						currentNode = triangleIndex;
					}
					else 
					{
						todo[todoOffset++] = bvhTree[currentNode].secondChildOffset;
						currentNode = currentNode + 1;
					}
				}
				//--
			} 
			else {
				if( todoOffset == 0 )
				{
					break;
				}
				currentNode = todo[--todoOffset];
			}
		}

		// Write result
		rayHits[gid].Distance   = maxT;
		rayHits[gid].U = b1;
		rayHits[gid].V = b2;
		rayHits[gid].Index = hitIndex;
	};
}